单细胞亚群鉴定知多少 第三篇
单细胞亚群鉴定知多少（三）
随着单细胞数据的增多，单细胞亚群的鉴定越来越成为一个单细胞分析的一个关键的、限制性的过程。
今天继续分享一下我们究竟如何去“读取单细胞这本天书”-单细胞亚群鉴定。
一
cellassign安装
cellassign（https://github.com/Irrationone/cellassign）是使用输入的表达矩阵，根据提供的marker基因，基于TensorFlow的概率模块进行机器学习的方法预测细胞类型。
因此安装的时候需要安装TensorFlow模块，也需要python，此工具最好用annoconda进行安装。
1 install.packages("tensorflow")
2 library(tensorflow)
3 install_tensorflow(extra_packages = "tensorflow-probability")
4 #Please ensure this installs version 2 of tensorflow. You can check this by calling
5 tensorflow::tf_config()
6 #> TensorFlow v2.1.0 (/usr/local/lib/python3.7/site-packages/tensorflow)
7 #> Python v3.7 (~/.virtualenvs/r-reticulate/bin/python)
8 sess = tf$Session()
9 hello <- tf$constant('Hello, TensorFlow!')
10 sess$run(hello)
11 BiocManager::install('cellassign'）
12 或者安装
13 devtools::install_github("Irrationone/cellassign")
14 ##
15 library(SingleCellExperiment)
16 library(Cellassign)
17 library(Seurat)
18 library(scran)
19 library(dplyr)
二
cellassign使用
cellassign数据的准备：
1 # 表达矩阵的准备
2 count<-as.matrix(rds@assays$RNA@counts)
3 # sizeFactors的准备
4 cell_anns <- data.frame(Barcode = names(Idents(rds)),celltype=Idents(rds),samples=rds@meta.data$stim)
5 rownames(cell_anns) <- names(Idents(rds))
6 sceset <- SingleCellExperiment(assays = list(counts = as.matrix(count)),colData=cell_anns)
7 str(sceset)
8 qclust <- quickCluster(sceset, min.size = 30)
9 sceset <- computeSumFactors(sceset, sizes = 15, clusters = qclust)
10 sceset <- normalize(sceset)
11 saveRDS(sceset,paste(oudir,paste(pref,'sceset.rds',sep="_"),sep="/"))
12 s <- sizeFactors(sceset)
13 # marker基因的准备
14 marker_mat <-read.table(markerfile,sep="\t",header=T)
15 marker_gene_list <- list()
16 for (i in marker_mat[,1]){
17 print(unlist(strsplit(as.character(marker_mat[marker_mat[,1]==i,2]),split = ",",fixed=T)))
18 marker_gene_list[[i]]<- unlist(strsplit(as.character(marker_mat[marker_mat[,1]==i,2]),split = ",",fixed=T))
19}
20 print(marker_list_to_mat(marker_gene_list))
21 marker_mat<-marker_list_to_mat(marker_gene_list)
22 fit <- cellassign(exprs_obj = sceset[as.vector(rownames(marker_mat)),],
23 marker_gene_info = marker_mat,
24 s = s,
25 learning_rate = 1e-2,
26 shrinkage = TRUE,
27 verbose = FALSE)
28
cellassign输入文件主要是三个，一个是表达矩阵，另外一个是marker基因list，最后一个是sizeFactors文件。
不过在这里sizeFactors计算直接通过rds文件的信息进行计算，因此没有单独列出，可以直接看上述命令即可（在进行computeSumFactors时候，有时候可能会报错，可以调大size参数即可）。
cellassign根据marker基因的list进行机器学习，预测细胞类型，因此此软件耗时较长，一般来说，10000个细胞的话，差不多24h。
三
结果展示
其结果展示也与之前讲的singleR一样，提供了热图展示，但是没有细胞亚群信息，因此我们需要把亚群信息添加上去，才更能反映我们数据的结果。
1#官方展示命令
2pheatmap::pheatmap(cellprobs(fit))
3
4cell_cluster_cellprobs<-t(cellprobs(fit))
5colnames(cell_cluster_cellprobs)<-sceset$Barcode
6annotation_col<-data.frame(samples=rds@meta.data$stim,Celltype=rds@meta.data$seurat_clusters)
7 rownames(annotation_col)<-rownames(rds@meta.data)
8 order_cells<-annotation_col %>% dplyr::mutate(.,barcod=rownames(.)) %>% dplyr::arrange(.,Celltype,samples)
9 gaps_col<-c()
10 m<-0
11 for (i in 1:max(as.numeric(annotation_col[,c('Celltype')]))){
12 gaps_col<-c(gaps_col,table(annotation_col[,c('Celltype')])[i]+m)
13 m<-table(annotation_col[,c('Celltype')])[i]+m
14 }
15 #热图展示命令
16pheatmap(cell_cluster_cellprobs[,as.vector(order_cells$barcod)],show_rownames=T,show_colnames=F,cluster_cols = F,cluster_rows= F,annotation_col = annotation_col,treeheight_row=0,border=FALSE,gaps_col = gaps_col,main = 'All Single-cell Recognition ',cellheight=16,fontsize=16)
cellassign最终会得到一个细胞类型的热图，颜色越深代表可能性越高。上图可以看出不同亚群、不同样品的细胞类型预测的可能性。颜色越红，预测可能性越高，最高为1。
http://www.bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html
一
clusterprofiler安装和使用
1if (!requireNamespace("BiocManager", quietly = TRUE))
2install.packages("BiocManager")
3BiocManager::install("clusterProfiler")
4
5cell_markers <- vroom::vroom('http://bio-bigdata.hrbmu.edu.cn/CellMarker/download/Human_cell_markers.txt') %>%
6 tidyr::unite("cellMarker", tissueType, cancerType, cellName, sep=", ") %>%
7 dplyr::select(cellMarker, geneID) %>%
8 dplyr::mutate(geneID = strsplit(geneID, ', '))
9cell_markers
10markerfile<-"all.markers.csv" #Seurat分析fFindAllMarkers函数得到所有亚群差异的输出csv文件
11seurat_marker<-function(markerfile,sep=",",colname="cluster"){
12 marker_list<-list()
13 print(paste("开始读取markerfile文件，时间为：",Sys.time(),sep=" "))
14 marker<-fread(markerfile,header=T,sep=sep,stringsAsFactors = FALSE,quote = "",data.table = F)
15 for (i in unique(marker[[colname]])){
16 marker_list[[as.character(i)]]<-as.vector(marker[marker[[colname]]==i,]$gene)
17 }
18 #names(marker_list)<-unique(marker[[colname]])
19 print(paste("完成markerfile转换成marker_list，时间为：",Sys.time(),sep=" "))
20 return (marker_list)
21}
22marker_list <-seurat_marker(markerfile)
23y <- compareCluster(marker_list, fun='enricher',TERM2GENE=cell_markers,minGSSize=1,pvalueCutoff = 1, pAdjustMethod = "BH", qvalueCutoff = 1)
24pdf(“compareCluster_cell_markers.pdf, w=12,h=8)
25p1<-dotplot(y, showCategory=as.numeric(showCategory),includeAll=TRUE)
26print(p1)
27dev.off()
28write.table(y@compareClusterResult,paste(outdir,paste(prefix,"compareCluster_cell_markers.xls",sep="_"),sep="/"),quote=F,sep="\t")
二
结果解读
这里可以使用enricher或者compareCluster函数进行富集分析，然后可以通过气泡图进行展示，最终得到如上所示气泡图，每一行为一个细胞类型，每一列为一个细胞亚群，颜色越红代表富集越显著性，气泡越大，代表富集程度越高。
目前来看，几乎所有的工具都是针对人和小鼠两个模式生物，较少的工具支持其他物种。
因此对于除了人、小鼠以外的物种而言，可能像scMCA/scHCA工具不能适用，marker基因数据库也有较多的限制，需要从文献中寻找想要的marker基因。
如果有参考数据集，可以根据已有的参考数据集通过singleR进行细胞类型预测或者根据marker基因通过cellassign和Garnett类似半监督的工具进行细胞亚群预测。
1. Abdelaal T, Michielsen L, Cats D, et al. A comparison of automatic cell identification methods for single-cell RNA sequencing data[J]. Genome Biology, 2019, 20(1): 1-19.
2. Han X, Wang R, Zhou Y, et al. Mapping the Mouse Cell Atlas by Microwell-Seq[J]. Cell, 2018, 172(5): 1307-1307.
3. Han, X. et al. Construction of a human cell landscape at singlecell level. Nature https://doi.org/10.1038/s41586-020-2157-4 (2020).
4. Park J, Shrestha R, Qiu C, et al. Single-cell transcriptomics of the mouse kidney reveals potential cellular targets of kidney disease[J]. Science, 2018, 360(6390): 758-763.
5. Zhang X, Lan Y, Xu J, et al. CellMarker: a manually curated resource of cell markers in human and mouse[J]. Nucleic Acids Research, 2019.
6. Franzen O, Gan L, Bjorkegren J, et al. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data[J]. Database, 2019.
7. Aran D, Looney A P, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage.[J]. Nature Immunology, 2019, 20(2): 163-172.
8. Zhang A W, Oflanagan C H, Chavez E, et al. Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling[J]. Nature Methods, 2019, 16(10): 1007-1015.
9. Yu G, Wang L G, Han Y, et al. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters[J]. Omics A Journal of Integrative Biology, 2012, 16(5): 284-287.
作者：尧小飞
审稿：童蒙
编辑：angelica
精选文章